#----------------------------------------------------------#
#  Democratic Resilience (Aurel Croissant and Lars Lott)   #
#----------------------------------------------------------#

# Title: "Democratic Resilience in the Twenty-First Century. Search for an Analytical Framework and Explorative Analysis" #
# Authors: "Croissant, Aurel and Lott, Lars", Heidelberg University and FAU Erlangen-Nürnberg
# date: 2025-05-07
# journal: Political Studies
# DOI: 10.1177/00323217251345779
# written under "R version 4.5.0 (2025-04-11 ucrt)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# load packages
library(here)
library(tidyverse)
library(magrittr)
library(reshape2)
library(ggpubr)
library(countrycode)
library(ggrepel)
library(ggthemes)
library(lubridate)
library(dotwhisker)

library(coda)
library(MASS)
library(rjags)
library(parallel)
library(magrittr)
library(dplyr)
library(tools)
library(devtools)

## Install local v-dem tools package
#devtools::install_local("vutils")
library(vutils)

# set you working directory here #

## load V-Dem data V14 ##

vdem_full_v14 <- readRDS(file.path("data/vdem14", "V-Dem-CY-Full+Others-v14.rds"))

## load Resilience Capacity data ##

ds.basic.subset.2000 <- readRDS("data/ds.basic.subset.2000_subset.rds")

## load ERT data v14 ##

ert_data <- readRDS("data/ert/ert_data_v14_uturns.rds")

summary(ert_data)

table(ert_data$uturn_ep_id)


bmr_data <- read.csv("data/bmr/bmr_data_extended.csv")

bmr_data <- bmr_data %>% dplyr::select(-X)

raw_data_indices <- readRDS("data/ds.basic.subset.2000.rds")

raw_data_indices <- raw_data_indices %>%
  dplyr::select(-c(country, iso2c, iso3c, iso3n,  historical_date, starts_with("v2x_polyarchy")))


#### Data Merging ####

vdem_subset_v14 <- vdem_full_v14 %>%
  dplyr::select(country_id, year, COWcode, country_name, starts_with("v2x_polyarchy"), starts_with("v2x_libdem"),
                starts_with("v2x_regime"), v2caassemb, v2cagenmob, v2caconmob, v2cademmob, v2caautmob, e_regionpol_6C,
                v2xca_academ, e_pt_coup, e_regionpol, v2x_regime, e_gdppc, e_pop, e_polity2, e_wbgi_pve)

ert_data <- ert_data %>%
  dplyr::select(-c(country_text_id, country_name, v2x_regime:v2x_polyarchy_codehigh, change_edi, first_v2x_polyarchy))

vdem_subset_v14 <- vdem_subset_v14 %>%
  left_join(ert_data, by =c("country_id", "year"))

vdem_subset_v14 <- vdem_subset_v14 %>%
  left_join(ds.basic.subset.2000, by =c("country_id", "year"))

vdem_subset_v14 <- vdem_subset_v14 %>%
  left_join(bmr_data, by = c("COWcode" = "ccode", "year"))

vdem_subset_v14 <- vdem_subset_v14 %>%
  left_join(raw_data_indices, by = c("country_text_id", "year"))

vdem_subset_v14 <- vdem_subset_v14 %>%
  distinct(country_id, year, .keep_all = T)

## Filter after 200 and only democracies ## 
vdem_subset_v14 <- vdem_subset_v14 %>%
  group_by(country_text_id) %>%
  filter(year >=2000) %>%
  filter(democracy == 1)

################################################################################
#### Internal Construct Validity Test ####
################################################################################

#### Macro-Institutional Dimension ####

cor_ResCap_v2x_horacc <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$v2x_horacc, 
                           use="complete.obs")
cor_ResCap_v2x_horacc  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_v2x_horacc), 2)))

## F1 ##
f1 <- ggplot(data = vdem_subset_v14, aes(x = v2x_horacc,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.2,
                 y = 0.9,label = cor_ResCap_v2x_horacc),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Horizontal Accountability Index",
       x = "Horizontal Accountability Index",
       y = "DRC Index")

## democracy Stock ##

cor_ResCap_democracy_stock <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$democracy_stock, 
                                  use="complete.obs")
cor_ResCap_democracy_stock  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_democracy_stock), 2)))

## F2 ##
f2 <- ggplot(data = vdem_subset_v14, aes(x = democracy_stock,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.2,
                 y = 0.9,label = cor_ResCap_democracy_stock),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Democracy Stock Index",
       x = "Democracy Stock Index",
       y = "DRC Index")

(f1 + f2) 

ggsave("../results/Figure_A6.png", height = 10, width = 15, units = "cm")


#### Political Actors Dimension ####

cor_ResCap_anti_pluralist_party_index <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$anti_pluralist_party_index, 
                             use="complete.obs")
cor_ResCap_anti_pluralist_party_index  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_anti_pluralist_party_index), 2)))

## F1 ##
f3 <- ggplot(data = vdem_subset_v14, aes(x = anti_pluralist_party_index,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.2,
                 y = 0.9,label = cor_ResCap_anti_pluralist_party_index),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Anti-Pluralist Parties Index",
       x = "Anti-Pluralist Parties Index",
       y = "DRC Index")


cor_ResCap_v2cacamps <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$v2cacamps, 
                                  use="complete.obs")
cor_ResCap_v2cacamps  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_v2cacamps), 2)))

## F2 ##
f4 <- ggplot(data = vdem_subset_v14, aes(x = v2cacamps,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.2,
                 y = 0.9,label = cor_ResCap_v2cacamps),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Political Polarization",
       x = "Political Polarization Indicator",
       y = "DRC Index")


cor_ResCap_v2caviol <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$v2caviol, 
                            use="complete.obs")
cor_ResCap_v2caviol  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_v2caviol), 2)))

## F2 ##
f5 <- ggplot(data = vdem_subset_v14, aes(x = v2caviol,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.2,
                 y = 0.9,label = cor_ResCap_v2caviol),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Political Violence",
       x = "Political Violence Indicator",
       y = "DRC Index")

(f3 + f4 + f5) 

ggsave("../results/Figure_A7.png", height = 10, width = 20, units = "cm")

#### Civil Society Dimension ####

cor_ResCap_v2xcs_ccsi <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$v2xcs_ccsi, 
                             use="complete.obs")
cor_ResCap_v2xcs_ccsi  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_v2xcs_ccsi), 2)))

## F1 ##
f6 <- ggplot(data = vdem_subset_v14, aes(x = v2xcs_ccsi,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.4,
                 y = 0.9,label = cor_ResCap_v2xcs_ccsi),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Core Civil Society Index",
       x = "Core Civil Society Index",
       y = "DRC Index")

## equal Access Index ##

cor_ResCap_v2xeg_eqaccess <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$v2xeg_eqaccess, 
                                  use="complete.obs")
cor_ResCap_v2xeg_eqaccess  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_v2xeg_eqaccess), 2)))

## F2 ##
f7 <- ggplot(data = vdem_subset_v14, aes(x = v2xeg_eqaccess,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.4,
                 y = 0.9,label = cor_ResCap_v2xeg_eqaccess),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Equal Access Index",
       x = "Equal Access Index",
       y = "DRC Index")

(f6 + f7) 

ggsave("../results/Figure_A8.png", height = 10, width = 15, units = "cm")


#### Political Community Dimension ####

cor_ResCap_Satis <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$Satis, 
                                             use="complete.obs")
cor_ResCap_Satis  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_Satis), 2)))

## F1 ##
f8 <- ggplot(data = vdem_subset_v14, aes(x = Satis,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.2,
                 y = 0.9,label = cor_ResCap_Satis),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Satisfaction with Democracy",
       x = "Satisfaction with Democracy",
       y = "DRC Index")


cor_ResCap_leg <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$leg, 
                            use="complete.obs")
cor_ResCap_leg  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_leg), 2)))

## F2 ##
f9 <- ggplot(data = vdem_subset_v14, aes(x = leg,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.2,
                 y = 0.9,label = cor_ResCap_leg),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Trust in Legislature",
       x = "Trust in Legislature",
       y = "DRC Index")


cor_ResCap_police <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$police, 
                           use="complete.obs")
cor_ResCap_police  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_police), 2)))

## F2 ##
f10 <- ggplot(data = vdem_subset_v14, aes(x = police,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.2,
                 y = 0.9,label = cor_ResCap_police),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Trust in Police",
       x = "Trust in Police",
       y = "DRC Index")

(f8 + f9 + f10) 

ggsave("../results/Figure_A9.png", height = 10, width = 20, units = "cm")

#################################################################################
#### External Validity Test ####
#################################################################################

## Liberal Democracy ##
cor_ResCap_v2x_polyarchy <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$v2x_polyarchy, 
                        use="complete.obs")
cor_ResCap_v2x_polyarchy  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_v2x_polyarchy), 2)))

## F1 ##
f11 <- ggplot(data = vdem_subset_v14, aes(x = v2x_polyarchy,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.3,
                 y = 0.9,label = cor_ResCap_v2x_polyarchy),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Electroal Democracy Index",
       x = "Electroal Democracy Index",
       y = "DRC Index")


cor_ResCap_v2x_libdem <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$v2x_libdem, 
                      use="complete.obs")
cor_ResCap_v2x_libdem  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_v2x_libdem), 2)))

## Electoral Democracy ##

f12 <- ggplot(data = vdem_subset_v14, aes(x = v2x_libdem,
                                         y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.3,
                 y = 0.9,label = cor_ResCap_v2x_libdem),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Liberal Democracy Index",
       x = "Liberal Democracy Index",
       y = "DRC Index")


(f11 + f12) 

ggsave("../results/Figure_A10.png", height = 10, width = 20, units = "cm")

## Polity2 ##
cor_ResCap_e_polity2 <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$e_polity2, 
                         use="complete.obs")
cor_ResCap_e_polity2  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_e_polity2), 2)))

## F2 ##
f13 <- ggplot(data = vdem_subset_v14, aes(x = e_polity2,
                                          y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = 0.2,
                 y = 0.9,label = cor_ResCap_e_polity2),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Polity IV",
       x = "Polity IV",
       y = "DRC Index")


cor_ResCap_e_wbgi_pve <- cor(vdem_subset_v14$ResCap, vdem_subset_v14$e_wbgi_pve, 
                            use="complete.obs")
cor_ResCap_e_wbgi_pve  <- paste("R =", sprintf("%.2f", round(mean(cor_ResCap_e_wbgi_pve), 2)))

## F2 ##
f14 <- ggplot(data = vdem_subset_v14, aes(x = e_wbgi_pve,
                                          y = ResCap)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = 'lm', formula = 'y ~ x') +
  geom_label(aes(x = -2,
                 y = 0.9,label = cor_ResCap_e_wbgi_pve),
             size = 5) +
  theme_clean() +
  ylim(0,1) +
  theme(legend.position="none",
        axis.text  = element_text(size=8),
        axis.title = element_text(size=9),
        plot.title = element_text(hjust = 0.5, size = 8),
        strip.background = element_blank()) +
  labs(title = "Political Stability (World Bank Governance Indicators)",
       x = "Political Stability",
       y = "DRC Index")


(f13 + f14) 

ggsave("../results/Figure_A11.png", height = 10, width = 20, units = "cm")
